The U(l) Gross- Neveu Model at Non-Zero Chemical Potential 



Simon Hands 

Department of Physics, University of Wales, Swansea, 
Singleton Park, Swansea, SA2 8PP, U.K. 



• Seyong Kim 

G^ 

Q>^ ' School of Physics, High Energy Physics Division 

^ ! 

^ ' Argonne National Laboratory 

9700 S. Cass Avenue, Argonne, IL 60439 

m 

John B. Kogut 

> : 

■ Department of Physics, University of Illinois at Urbana- Champaign 

o 



o 



- 

X 



1110 West Creen Street, Urbana, IL 61801-3080 



Abstract 



■ The four-fermi model with continuous chiral symmetry is studied in three di- 



mensions at non-zero chemical potential /i using both the ^/Nf expansion and 



■ computer simulations. For strong coupling this model spontaneously breaks 

its U(l) chiral symmetry at zero chemical potential and the Goldstone mech- 
anism is realized through massless pions. The computer simulation predicts a 
critical chemical potential close to the lightest fermion mass in the model. 
As )U is increased beyond the pion screening mass increases rapidly from 
zero to a nonvanishing value indicating symmetry restoration. Some lessons 

are drawn relevant to lattice QCD simulations at non-zero /x. 
ll.lO.Kk, 11.15.Ha, 11.15.Pg, ll.lO.Wx 
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I. INTRODUCTION 



Simulations of lattice QCD at non-zero chemical potential /i have remained in a quagmire 
for over a decade The lattice action becomes complex when /i is taken non-zero, so 
conventional computer simulation algorithms which are based on a probability distribution 
do not apply. The quenched version of lattice QCD, which attempts to skirt the issue by 
replacing the complex fermion determinant by unity, appears to be pathological in the chiral 
limit [^]. When the bare quark mass is non-zero so the pion is not massless, one finds 
that quenched simulations of lattice QCD work for /i outside a "forbidden" region extending 
from to m^/S, where m^^^ denote pion and baryon masses respectively The failure 

of the quenched version of QCD to describe the forbidden region is poorly understood ||5| 
[0. Studies of zero dimensional models of non-zero chemical potential have not been 
decisive - some models work in their quenched versions [§| and others do not 0. Since such 
models cannot respect the Goldstone mechanism, their relevance to QCD is questionable at 
best. 

In this paper, we study the A^j-fiavor four-fermi model with U(l) chiral symmetry, some- 



times called the Gross- Neveu model [|T0|, in three dimensions (ie. two space and one time) 
as a function of chemical potential fi. The action is 



c = i>i{^ + /i7o + m)i)i - ^^[{AAY - {iPa5Ay], (i-i) 

where ■?/','?/' are four component spinors and a sum on i = 1,...,A^/ is understood. For 



fi < jjc this model has a massless pion in the chiral limit m — ||rT| and can be studied by 
the leading order expansion [l^ [T^ as well as by conventional computer simulations. 
Since the pathologies in QCD simulations are tied to m^r, we can address the possibility that 
the presence of this massless mode is the culprit in past failures. In addition, the model 
has composite mesons like QCD and has sensible infra-red and ultra-violet properties. The 
zero chemical potential model has a chiral transition as a function of its bare coupling which 



can be analyzed using the expansion [|T5[. The second order phase transition 
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coincides with a renormalization group fixed point and defines an interacting continuum 



field theory [|I2| |T^. It will prove interesting to simulate the U(l) four-fermi model to see 
the physics of chiral symmetry restoration at work in a model with a realistic, composite 
pion. The four-fermi model lacks two central features of QCD: 1. the model does not 
confine, and 2. its fermion determinant is real and non-negative even when /i 7^ 0. 

We shall see that the simulation study of the four-fermi model is completely successful 
and physical, and in agreement with the predictions of the expansion. In addition to 
the usual thermodynamic quantities such as the chiral condensate (ipip), the fermion energy 
density e, and the fermion density p itself, we will measure some spectroscopic features 
of the model. These will include the pion and fermion screening lengths, since they are so 
closely tied to the physics issues of interest. Recall how chemical potentials are implemented 
and screening lengths are measured. Consider a symmetric lattice and label one of the 
axes r, "temporal." In the +r direction assign a factor exp(/i) to each such link in the 
action, and in the — r direction assign a factor exp(— /i). These exponential factors favor 
quark propagation in the +r direction and when fi is small in units of the reciprocal lattice 
spacing, ie. fia 0, we have a useful, well-behaved transcription of the chemical potential 



to a discretized system |T^ |[T^. Screening lengths are calculated in this environment by 
calculating propagation with a source at = and a sink at tj = r. The exponential fall- 
off with r then gives an estimate of the appropriate screening length, as will be illustrated 
through detailed calculations in the text below. The simplest expectation is that the pion 
inverse screening length should be zero in the chiral limit for small fi where chiral symmetry is 
spontaneously broken, and as /i increases through a critical value /ic, where chiral symmetry 
is restored, the pion inverse screening length should increase from zero non-analytically. 
The fermion inverse screening length should decrease linearly with fi and vanish aX = 
where is the value of the fermion mass at /i = 0. Interactions are expected to decrease /ic 
below this free field result, but at large Nf the naive result should be accurate. The result 
fic = assumes that the fundamental fermion is the lightest fermion in the model's mass 
spectrum. Although this result is expected in three dimensions, it is not true in four-Fermi 



models in two dimensions, where there are "kink" solutions and kink - anti-kink bound 



states 18 . Computer simulation of the two dimensional model were successful in finding 



the kink solutions and the subtle theoretical predictions were confirmed |T9[. In this study 
we shall confirm the naive expectation fic = rrif with good control. In fact, the /^-dependence 
of the fermion screening length will be consistent with, 

m/(/i) = m/(0) -/i (/X < /ic = m/(0)) (1.2) 

which is also the free-field prediction. 

The results of our 1/A^/ calculation and lattice simulations were very clear. The simula- 
tion gave the expected physical answers for both large Nf {Nj = 12), and intermediate Nf 
{Nf = 4). We worked near the chiral limit and found no pathologies when /i passed through 
the value m^/2. In fact, the algorithm was well-behaved for all fi and computer simula- 
tions in the immediate vicinity of fic = ""^/(O) were only limited by finite size effects of the 
expected variety. The chiral symmetry restoring transition was clear in all the observables 
calculated and excellent consistency was found between the various measurements. 

Does this success help us understand the confused state of QCD simulations at non-zero 
chemical potentials? It certainly indicates that a massless pion is not the culprit behind the 
failures encountered in QCD simulations. The compositeness of the pion is also seen to be 
harmless. As suggested by many authors in the past the complex nature of the QCD 
fermion determinant must be the source of the problem. The quenched version of QCD 
ignores the complex part and possible rapid variation in the determinant and is probably 
pathological because of this, as illustrated in the single-site U(l) model of Gibbs 0. A 
successful simulation of QCD at non-zero fi may require a wholly new algorithmic approach. 
The present generation of fermion algorithms calculate the fermion propagator in a given 
gauge field background, and this intermediate step may not yield useful results in QCD for 
m^/2 < /i < uib/S g. 

This paper is organized as follows. In Sec. II we define the lattice version of the U(l) 
four-fermi model and present some leading order large- iV/ results for the first order chiral 



symmetry restoring transition for // > in the model, which will be compared to computer 
simulations later in the paper. In Sec. Ill we consider the computer simulations of the 
Nf = 12 and 4 models. We shall find that the usual bulk thermodynamic observables such 
as (ipip) , energy densities and number densities successfully expose the chiral transition. In 
addition, spectroscopic measurements are equally clear and particularly informative. The 
pion inverse screening length is essentially zero in the broken, low phase, but is non-zero 
in the unbroken phase. The fermion inverse screening length decreases linearly with /i in 
the broken phase and vanishes at the same critical point seen in the other observables, 
lie — mf{0), the fermion mass observed at = 0. The hybrid Monte Carlo algorithm we 
use is efficient for all /i and the convergence of its underlying conjugate-gradient routines to 
invert the lattice Dirac operator does not deteriorate with increasing fj,. Sec. IV includes a 
brief summary and remarks on the lessons learned in this study. 



II. LATTICE FORMULATION OF THE GROSS-NEVEU MODEL 

The lattice action for the bosonized Gross-Neveu model with U(l) chiral symmetry is 

1=1 '-x,y ° X <x,x> <x,x> 

+ ^E(^'(^) + ^'(^-))- (2.1) 

Here, Xi and Xi are complex Grassmann- valued staggered fermion fields defined on the lattice 
sites, the auxiliary scalar and pseudoscalar fields a and tt are defined on the dual lattice 
sites, and the symbol < x,x > denotes the set of 8 dual sites x adjacent to the direct lattice 
site X. The auxiliary fields only appear to quadratic order, and can be integrated out to 
recover a lattice action in terms of fermion fields only (bosonization is also possible for the 
continuum Lagrangian (1.1)). The parameter Nf will turn out to be the number of physical 
fermion species described by the numerical simulation, and must be an integer multiple of 
4. The fermion kinetic operator Ai is given by 
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^y,x+u ^y,x—u 



+ m6. 



y,xi 



(2.2) 



where m is the bare fermion mass, /i is the chemical potential, and riy{x) are the Kawamoto- 
Smit phases (— l)'^'""' t"^"-!. The symbol e{x) denotes the alternating phase (_ l)^o+a:i+a;2_ 

As described in refs. [jl^ |T^, the model (2.1) can be recast in terms of fields q{Y) defined 
on a lattice of twice the spacing, where the field q has explicit spinor and flavor indices. In 
momentum space the kinetic part of the action reads 



S, 



kin 



E / 7^ E -fe(/c)(7.®l2)gi(A;)sin2A;, + gi(A;)(74®r;)g,(A;)(l-cos2A;,)] 



1 

+ 2 



qi{k){-yo (8> l2)qi{k)[i sin2/co cosh fi + {1 + cos 2/co) sinh fj] 



(2.3) 



+Q'i(^)(74 ® ''"3)Q'i(^)[^(l ~ cos2A;o) cosh/i + sin2A;o sinh/i] 



+ mqi{k){li ®l2) qiik), 



with Tj the Pauli matrices, which act on the 2 component flavor degrees of freedom, and 

• (2.4) 



7u 



7o 



74 



75 



-T-^/ V y \^l2 / VI2 

The momentum integral extends over the range k G (— 7r/2, 7r/2]. On a flnite system the 
integral is replaced by a sum over L/2 modes, where L is the number of lattice spacings in 
a given direction. In the large wavelength limit A; ^ we recover the standard Euclidean 
continuum form qj{<^ + fi'jQ + m)qj, where j now runs from 1 to Nf/2. 
The interaction part of the action can be rewritten as 



Sint = E HyMy)i^4 ® l2)q^iY) + 7r(f)g,(F)(275 ® l2)qiiY) + 0(c 



(2.5) 



Y 



where d'{Y), tt(Y) denote the sum of the eight scalar flelds associated with the site Y on 
the blocked lattice (see fl^ for full details), and a is the lattice spacing, which has been set 
to unity in Eqns. (2.1-3). The 0(a) terms are non-covariant and flavor non-singlet - if we 
used a formulation in which the scalar flelds were deflned on the direct lattice sites, then 
such non-covariant terms would appear at 0{a^) [pC|] . 

The lattice action (2.1) has a global vector U(A^//4) symmetry 



Xi 



(2.6) 



and, in the chiral limit m — 0, a global "axial" U(l) symmetry 

X ^ e'^<^\ ; X ^ Xe^"^^"^ ; cj^ ^ {a + in) ^ e'^^^. (2.7) 

In the g-basis the rotation (2.7) reads: 

g 1-^ expza(75 ® l2)g ; g i-^ gexp ^(75 (g) I2). (2.8) 

The U(l) symmetry is spontaneously broken by the condensate (xx) (or equivalently (cr)): 
for m 7^ this direction of symmetry breaking is picked out. 



The action (2.1) may be simulated using a hybrid Monte Carlo algorithm ||2T|, in which 
the Grassmann fields are replaced by complex commuting pseudo-fermion fields il){x) gov- 
erned by the action 

S = T.ll :^V^r(^)(MtM);i^.^,(2/) + -4EK(5) + ^'(5^))' (2-9) 

xy ij=l ^ °9 X ^ ^ 

where 

Mxyi2 = MxySij + SxySij^ [a{x) + i€{x)n{x)]. (2.10) 

<x,x> 

Integration over the ip fields yields the functional measure det(M'''M). Note that the kinetic 
part of M, Ai, is strictly real even for /i 7^ 0, and that the complex part a + isTC is diagonal. 
Thus, schematically, 

det (M^M) = detMMetM 

= det{M + a + ien) det{M + a - ien). (2-11) 

Since each matrix Ai effectively describes Nf/2 fermion species, we conclude that the hybrid 
Monte Carlo simulation describes a system of Nf fermion species, Nf/2 having positive chiral 
charge and Nf/2 negative. To be more precise, the full symmetry of the lattice model in the 
continuum limit is U(A^//2)v®U(iV//2)y (g)U(l)A rather than the U(A^/)y (8)U(1)a we would 
obtain if all species transformed identically with respect to chiral rotations. At non-zero 
lattice spacing, of course, the symmetry group is smaller as discussed above: U(A^//4)v ® 



\J{Nf/4)v ® U(1)a- In all cases, of course, it is the U(l)^ symmetry which is broken, either 
spontaneously by the dynamics of the system, or explicitly by a bare fermion mass. 

It was found that the performance of the hybrid Monte Carlo procedure could be op- 
timized in two ways: firstly, as described in [|T^ |T^, by tuning the effective number of 
fermion species N'jr used in the guidance part of the program (ie. during the integration of 
the equations of motion along a microcanonical trajectory) to maximize the acceptance rate 
of the Monte Carlo procedure for fixed timestep At; and secondly by choosing the trajectory 
length T at random from a Poisson distribution of mean f. This second refinement, which 
guarantees ergodicity, was found to dramatically decrease autocorrelation times p2| . 

In our specific application we found that the choice At = 0.01 was adequate - the 
acceptance rate in the Monte Carlo procedure was always high (better than 80%) while the 
algorithm sampled configuration space with good speed and efficiency. For the Nf = 4 model 
the acceptance rate was improved by taking Nj^ slightly larger, typically 4.05, although the 
"best" choice of N'j/Nf certainly depends on lattice sizes (16'^ in our case) and couplings 
(0 ^ 1/(7^ ^ 1.0) and bare fermion masses (0.01 in lattice units). The trajectory lengths 
were chosen from a Poisson distribution with r typically between 1 and 2. In the immediate 
vicinity of the critical /i, a larger r probably would have been better, but we had no trouble 
obtaining good data with modest error bars from runs of several hundred trajectories. The 
data and error bars will be presented in tables in the next section. We used identical 
parameters when the Nf = 12 model was simulated except N'j was now tuned typically to 
12.12, for better acceptance rates. 

As well as measuring the expectation value of the scalar field (a) in the simulation, which 
for our purposes serves as an order parameter, we also monitored the chiral condensate {xx)i 
the energy density (e), and the fermion number density (p), which are defined as follows: 

- {XX) = \^^^Sf = ^{irM-^), 
(e) = = LtTdoioSp = —iYe^M-^ , - e-^M"^ (2.12) 
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Here Vg is the spatial (ie. two-dimensional) volume, P the inverse temperature, and V = K/3 
the overall volume of spacetime. The final expression in each case is the quantity measured 
in the simulation, using a noisy estmator to calculate the matrix inverses. 

Finally, as well as results from numerical simulations, we can examine the action (2.1) 
using the expansion. To leading order, this corresponds to the saddle point solution 

of the path integral, or equivalently to mean field theory. We can solve for the expectation 
value of the scalar field (cr) using the gap equation [0 [|1^]: 

{a) = -g\xx) = ytrSF. (2.13) 

To leading order the gap equation consists of a single tadpole. Using the form for the kinetic 
term (2.3), we solve self-consistently to get an equation for in terms of the dynamical 
fermion mass fnf = (a) + m. In the thermodynamic limit L — > oo, 
1 _ 8 rrif r^n d?k 



g"^ 2ti^ rrif — m J -it/2 | { 1 — cos 2ko cosh 2fj, — i sin 2ko sinh 2fj,} + J2iy=i,2 sin^ + mj- ' 

(2.14) 

Note that the 0(a) terms in the interaction (2.5) make no contribution at this order. Eqn. 
(2.14) can be reduced to a one- dimensional integral and then evaluated numerically in the 



limit /i ^ [Q: however for /i 7^ it is more efficient to simply evaluate the sum over 
lattice momenta on finite systems explicitly. Note that antiperiodic boundary conditions 
must be specified on the fermion fields in the timelike direction; we chose to apply periodic 
boundary conditions in the two spacelike ones. 

In Figs. 1 - 5 we show the large- A^j predictions from Eqn. (2.14) evaluated on a 16^ 
lattice. In Fig. 1 we show (cr) vs. l/g"^. We have chosen bare fermion masses of m = 
0.05, 0.01 and 0.00. This range of masses shows the sensitivity of the finite lattice results to 
m, gives the theoretically interesting chiral limit m = 0.00 and includes m = 0.01 to compare 
quantitatively to the simulation results obtained from the hybrid Monto Carlo algorithm at 



Nf = 12 and 4,m = 0.01. A chiral symmetry restoring symmetry point at 1/gi slightly 



larger than 1.00 is found consistent with the critical index (3.mag = 1-0 The curves at 
m = 0.01 and 0.05 smooth the transition out in the expected fashion. Only the m = 0.00 
curve shows evidence of finite volume rounding (this was checked by evaluating the gap 
equation on larger lattices). Note that (cr) is also the fermion dynamical mass mj in the 
chiral limit, so mj is an equally good order parameter for the transition. Since mj is an 
inverse correlation length, the critical exponent u = 1.0 follows from the figure. In Fig. 2 we 
show (cr) vs. fi for a bare coupling 1/ g"^ = 0.5 deep within the broken symmetry phase. This 
coupling will be simulated in the next section. We note that (a) is essentially unaffected by 
11 until the immediate vicinity of the transition where (cr) jumps to zero through a first order 
phase transition. (Recall that in mean field theory a model can have a phase transition in a 
finite volume.) The rounding of the curves near /x = /i^ is a finite volume effect. The critical 
/ic is approximately 0.73 which compares rather well with the dynamical fermion mass in the 
chiral limit at l/g"^ = 0.5 and fi = 0.0, = 0.84, seen in Fig. 1. The discrepancy between 
He and is simply due to the discreteness of the lattice sums approximating continuum 
integrals. In the continuum limit g ^ gc (ie. ^ 0), the discrepancy would approach zero. 
It is interesting to note that if the plot of mj vs. l/g"^ in Fig. 1 is extrapolated back from 
1/g'^ using the assumption of linearity (ie. Pmag = 1), then the resulting rrif at l/g"^ = 0.5 
is in much closer agreement. The m = 0.01 and 0.05 curves in Fig. 2 demonstrate that 
these bare masses do not distort the first order symmetry restoring transition significantly. 
Since the mass ratios m/mf are 0.012 and 0.058, approximately, in the two cases and since 
the transition is strongly first order, this result is expected. In Fig. 3 we show the lattice 
gap equation prediction for a vs. /i at m = 0.01 on a 16'^ lattice for values of l/g"^ ranging 
from 0.5 to 0.8, showing the effect of finite volume rounding becoming more pronounced 
as the physical lattice spacing is reduced. By 1/g'^ = 0.8 the transition has become much 
less sharp. In Figs. 4 and 5 we show the fermion number density p and the fermion energy 
density e plotted against fi in the chiral limit for three choices of couplings in the broken 
phase, l/g"^ = 0.5,0.7 and 0.9. Also shown in Fig. 4 is the expected continuum result for 
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fermions in the unbroken phase, p = p^/n (we have rescaled the numerical results for p 
and e by a factor of two over the definitions (2.12)) [|1^]. For l/g"^ = 0.9 the transition is 
hard to identify. As we approach the bulk critical point 1/ g1 1.0 on the finite system the 
signals of the phase transition become less and less dramatic. Larger lattices are necessary 
to obtain a quantitative picture of the transition when 1/ g"^ is chosen to be 0.9, close to the 
bulk continuum limit 1/ g1 ~ 1.0. 



III. SIMULATION RESULTS 

We simulated both the Nf = 12 and Nf = 4 models at various couplings and chemical 
potentials. In all 16^ lattice was used, with a bare fermion mass m = 0.01. As 

explained above, due to fermion "doubling," the Nj = 12(4) model corresponds to A^/iatt = 
3(1) lattice species. We simulated the A^/iatt = 3 model because it should compare well 
with the results of the l/Nj expansion (if its underlying assumptions hold in the model), 
and observables should be large with relatively modest fluctuations. The A^/iatt = 1 model 
is interesting because it could show qualitative deviations from the laige-Nf results; its 
observables will fluctuate more intensely and it will present a numerical challenge closer 
to that of two- or four-flavor QCD. The observables and measurement techniques were 
discussed above. They have also been used in our past studies of four-fermi models, so 



we refer the reader to those references rather than repeat standard material [|1J]. The 
"new" measurements we have done concern the model's spectroscopy. These observables are 
particularly revealing and we will discuss them at greater length. 

Using point sources, we calculated the fermion propagator, G^{x,t), with chemical po- 
tential /i and, G-{x,t), with chemical potential —p. Then we formed a zero momentum 
fermion propagator 

Pf{t) = T.G+{x,t), (3.1) 

X 

and anti-fermion propagator 
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Pfit) = J2G4x,t). (3.2) 

X 

Of course, t) = G^{x, t) for the zero chemical potential case. The composite pion and 

sigma propagators are, 

PAt) = T.G+iS,t)GUx,t), (3.3) 

X 

P^t) = J2{-ir+yG4x,t)Gl{x,t), (3.4) 



in analogy with the treatment using staggered fermions in four dimensions |2^. We also 
calculated propagators for the auxiliary fields in Eq. (2.1), vr and a, 

^7r(t) = 5I^(^.^')7r(x,t' + t), (3.5) 

x,t' 

and 

P^{t) = 5][(T(f, t')a{x, t' + t)- a\ (3.6) 

x,t' 

Here, is the square of the average of the a field. 

After calculating averages of the above propagators and their covariance matrices (see 
p3| for fitting techniques), we fit the various propagators to the following functional forms. 
For zero chemical potential: 

Pf(t) = A[e-™/* - (-l)*e-"^^(^-*)], (3.7) 

where T is the temporal extent of the lattice, for the fermion and 

p^{t) = Ale-""-' + e^™-(^-*)], (3.8) 

for the auxiliary field, vr. Similar functional forms are chosen for the anti-fermion and a 
respectively. For the non-zero chemical potential case, we chose the forms, 

Pf{t) = A[e-'^'f' - (-l)*e-'"/(^-*)] (3.9) 

P^it) = A[e-"^^* + e-'"'(^-*)]. (3.10) 
12 



Since the meson has zero fermion number, we expect that ml — m^. 

Those parameters which minimize correlated were chosen as the best fitted values. 
We used the CERN mathematical library routine, MINUIT, as a minimization program. 
The error bars quoted refer to the necessary parameter changes for a change of by one. 
Only two of our spectrum calculations yielded useful, accurate mass estimates. They were 
the fermion and the auxiliary field pion masses. Luckily, these are the two quantities most 
closely related to chiral symmetry and its restoration at non-zero chemical potential. It 
would have been interesting to calculate the a mass, but the fluctuations in our limited data 
sets made that impossible. 

Now let's turn to the data. The A^/iatt — 3 data for m/,m7r, the vacuum expectation 
value (a) , which we shall denote a for simplicity in this section, and the action S are given 
in Table I for /j, = and couplings ranging from 0.50 to 1.0. The order parameter a and 
ruf are almost identical as they should be at large Nf. We show a plot of a vs. l/g"^ in Fig. 
6. Clearly chiral symmetry is broken for ^ 0.90 — 1.00. Also shown is the prediction 
of the lattice gap equation from Fig. 1. The agreement is 0(10%) or better, and is very 
satisfactory. It is our first indication that the 1/A^/ expansion is practical and obtains the 
correct physics of these models. We note from the table that the mass of the pion is "small" 
and quite insensitive to 1/ g'^. Since chiral symmetry is broken over this range of couplings, 
it should be that the pion's mass is nonzero only because of the explicit symmetry breaking 
provided by the small bare fermion mass, m — 0.01. We will obtain good evidence for this 
interpretation of the data when we consider the model at nonzero chemical potential. 

We studied the model next at l/g"^ — 0.50 as a function of // to see how a nonzero 
chemical potential restores chiral symmetry at a critical point. \jg^ — 0.50 is a good place 
to look first because it is deep in the broken symmetry phase, but not too far from the bulk 
critical point l/g"^ 1.00. As long as l/g"^ ^ 0.50 is in the scaling window of the critical 
point we can extract continuum physics from these simulations by standard methods. Our 
goals here are more modest - we wish to turn up n and see if the conventional picture of 
symmetry restoration emerges. The simulation data for 1/g^ — 0.50 and /i ranging from 
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0.50 to 0.85 is given in Table II. The table includes data for ruf (the fermion mass), m.,^ 
(the pion mass), p (the fermion number density), a (the vacuum expectation value of the 
a field), S (the action), and e (the energy density). Since the fermion mass at /i = 0.0 is 
rrif = 0.746(2) according to Table I, we expect naively a chiral symmetry restoring transition 
at /ic = 0.746(2). This result is beautifully reproduced by the simulation. In Fig. 7 we plot 
the order parameter a vs. fi for fixed l/g"^ = 0.50, and we see restoration of the symmetry 
at pc = 0.725(25). The curve is very abrupt and we suspect, naturally, that a simulation 
on a larger lattice at smaller bare fermion mass m would show a first order discontinuous 
transition. The curve is also quantitatively very similar to the predictions of the lattice gap 
equation which is also shown - once again, the discrepancy is 0(10%) and may presumably 
be ascribed to 0{1/Nf) corrections. 

As discussed in the introduction, one of the purposes of studying this model was to verify 
that the lattice formulation, given a proven algorithm, obtains the correct physics at non- 
zero chemical potential even in a model with a Goldstone pion. Quenched QCD simulations 
have pathologies when fi approaches m^/2, and although the chiral restoring transition is 
expected at fic = "^b/S, one-third the mass of the nucleon, there is little numerical evidence 
for this [0 01 [Q- Many reasons have been proposed in hindsight for this catastrophe and 
several of them hinge on the presence of a Goldstone pion in the theory's spectrum. In the 
simulations here 171^^/2 0.09 while the expected transition is at pc = = 0.74(2). In 
QCD simulations the two mass scales, and nis/'i-, are actually quite close in numerical 

simulations performed to date, and this has further clouded the situation. In our U(l) four- 
fermi model, we can simulate the model very close to the chiral limit in the sense that the 
explicit breaking is much smaller than the dynamical breaking (ie. bare mass m <^ rrif, the 
dynamical mass). Fig. 7 and the rest of the data in Table II show that 171^^/2 is not a point of 
any special significance and the simulation algorithm gives the expected pc = 0.74(2) nicely. 
Additional simulations show that the value m = 0.01 is not crucial to these conclusions, and 
the chiral limit m ^ does not contain surprises even when fi ^ 0. For example, in Fig. 8 
we show a vs. m for l/g"^ = 0.50 and fi = 0.50. Unlike the tortuous situation in quenched 
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QCD where plots of {ipip) vs. m have strong downward curvature for /i ^ "^tt/S [||], the 
extrapolation to the chiral limit here is essentially linear and without surprises. Simulations 
at larger values of m would have been as clear as the m = 0.01 case studied here in detail. 

The transition at /ic = 0.74(2) is seen equally well in the fermion number and energy 
densities (Fig. 9), and the action itself (Fig. 10). These curves compare well with the 
analytic large- A^/ results presented in Sec. II above. The spectroscopy of the model and its 
dependence on /i is particularly interesting. A recent study of quenched QCD |^ showed that 
in the limited range < fi < the baryon mass as defined through the exponential falloff 

of a Euclidean propagator as in Eq. (3.7) decreased linearly with fi and would vanish by linear 
extrapolation at /ic = ms/^ as expected. Unfortunately, the quenched simulation algorithm 
suffers from slow convergence and large fluctuations for ^ ^ w-tt/S, so nothing is known 
quantitatively in the "forbidden" region, < fi < tub/ 3. In addition, the quenched 

QCD simulation showed that the pion mass, as defined through the exponential fall-off of 
a propagator as in Eq. (3.8), is insensitive to fi for fi ^ This is another sensible 

result which could not be confirmed at larger /i due to the pathologies of the quenched 
simulation. In our four-fermi model the analogous calculations are successful for all fi. 
In Fig. 11 we show the fermion mass as obtained from Eq. (3.7). Up to modest and 
expected finite size effects which reduce the fermion mass estimates in the vicinity of fic, 
the calculation is successful, and gives a critical chemical potential near 0.74(2), although 
larger lattice studies and a systematic analysis of finite size effects would be necessary to 
obtain a quantitative prediction. Also in Fig. 11 we show the pion mass m^r, calculated 
using Eq. (3.8), as a function of fi. The pion mass exhibits no fi dependence until we reach 
the vicinity of fic = 0.74(2) where it jumps up indicating chiral symmetry restoration. We 
wanted to measure the a mass over the same range, but its propagator was noisier than 
the pions and quantitative estimates were not achieved. We had hoped to verify that the 
pion and the sigma are degenerate and heavy for fi > 0.74(2), indicating chiral symmetry 
restoration. Although from this perspective our calculations were only partially successful. 
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they gave decisive physical answers expected of a Goldstone particle as we pass through a 
chiral symmetry restoring transition. 

Next we repeated these measurements at l/g"^ = 0.80 in the A^/iatt = 3 model in order 
to be closer to the model's continuum limit. At l/g"^ = 0.80 the fi = value of the fermion 
mass is rrif = 0.295(3). On the basis of the analytic large- A^/ results shown in Figs. 3,4 
and 5, we anticipate that the transition will be harder to identify since finite size effects are 
more severe and smooth the transition considerably on a 16^ lattice. We show the data for 
= 0.80 in Table III. It is organized just as Table II was. In Figs. 12 - 14 we show 
0", p, e and S plotted vs. /i. These figures show the same features as the large- A^/ calculation, 
and indicate that the transition is in the vicinity of /ic = 0.295(3). There is no evidence 
whatsoever for pathological behavior at /i = m^/2 0.11(1). The transition is shown with 
greater clarity in the model's spectroscopy. In Fig. 15 we see that the dynamical fermion 
mass decreases, up to the expected finite size effects, linearly with p and vanishes when /x 
becomes 0.25(5). Similarly, the pion mass is insensitive to p until p = 0.30(2), where it 
increases noticeably. We learn that the model's spectroscopy is a more sensitive guide to 
the chiral restoration transition than traditional bulk thermodynamic quantities. In light of 
our recent work on quenched QCD , this result comes as no surprise, but it reinforces the 
strategy we are taking in the four dimensional gauge model. 

Next we turn to the A^/iatt = 1 model to simulate a case where fluctuations are expected 
to be more significant, the number of flavors is more realistic and the large- A^j expansion 
may not be as good a guide. The p = data is collected in Table IV. We notice, as 
plotted in Fig. 16, that the order parameter a and dynamical fermion mass nif disagree 
slightly deep in the broken symmetry phase, but are otherwise in good agreement. Since 
these two quantities are identical at large- A^j, we see signs of 1/A^/ corrections here, but 
they are not numerically significant near the transition. We will investigate the theory at 
non-zero yU at both \jg^ = 0.60 and 0.70. The Xjg^ = 0.60 simulation is quite far in the 
broken phase and should be relatively decisive while the Xj = 0.70 simulation will be more 
strongly affected by fluctuations. In Figs. 17 - 19 we show o", p, e and S plotted against \x 
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using the data of Table V. In this case the plots show only qualitative agreement with the 
predictions of the expansion at leading order; 0(1/A^/) corrections are numerically 

much more significant. We expect a transition near the value of the dynamical fermion 
mass at — 0.60, ie, m/ = 0.475(5), and the figures are in fine agreement with that. 
In addition the fermion and pion masses. Fig. 20, show the transition almost as clearly 
and quantitatively as they did in the A^/iatt = 3, = 0.50 case. Finally in Table VI 
and Figs. 21-24 we show the analogous quantities for the A^/iatt = 1 theory at l/g"^ — 0.70. 
Although the bulk thermodynamic quantities experience considerable rounding, the results 
are consistent with a critical chemical potential Hc = 0.32(1), as predicted by the value of 
the dynamical fermion mass at — 0.70. Once again, the spectroscopic quantities in Fig. 
24 provide more quantitative information. All in all, the simulations are successful in each 
case and do not suffer from the pathologies affecting quenched QCD. 



IV. CONCLUDING REMARKS 

The success and clarity of these simulations shows that the Hybrid Monte Carlo algorithm 
(and hence presumably the closely related Hybrid Molecular Dynamics algorithm, suitable 
for values of Nf which are not a multiple of 4) are completely rehable for this class of fermion 
field theories in which there is a massless pion in the chiral limit. A chiral symmetry restoring 
phase transition, probably first order, is found for a critical value of the chemical potential ii. 
Screening length calculations proved to be particularly illuminating and the expected physics 
of chiral symmetry restoration at the critical chemical potential emerged. The spectroscopic 
data may well prove to be the most accurate means of determining the critical chemical 
potential on finite systems. The predictions of the large Nf expansion proved to be a good 
guide into the physics of these four fermi models even when Nf assumed modest values. 
Systematic effects are of the expected form, and we have no reason to suspect they could 
not be brought under complete control given sufficient computer time. 

One of our primary motivations for this work was to narrow down the source of difficulties 
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in simulations of lattice QCD. We certainly have demonstrated that the presence of a mass- 
less pion in the theory's spectrum is not the source of those difficulties. By default it must 
be the complex nature of the QCD action at nonzero fi that is the culprit. The quenched 
version of QCD ignores the fermion determinant and this omission apparently amounts to 
a qualitative error when the chemical potential is larger than half the pion mass . In four 
fermi models the fermion determinant is real and non-negative and the Euclidean theory 
has a probabilistic interpretation. In addition, the fermion propagator is well-behaved for 
all /i. These ingredients then allow the Hybrid algorithms to be successful in four fermi 



models despite being inapplicable for QCD. It may well be that Langevin algorithms |24 



can simulate QCD directly at nonzero /x, but that is far from clear at this time (Langevin 
algorithms are known to correctly evaluate certain complex integrals, and certain physical 
systems with complex actions, but their successes and failures are difficult to anticipate in 
advance). It will probably require greater insight into the numerics of the complex fermion 
determinant before a truly trustworthy, first principles simulation of QCD in an environment 
rich in baryons is possible. 
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TABLES 

TABLE I. Ajiatt = 3 data on 16^ lattice at /x = 0.0. The columns give the coupling 1/^^, the 
fermion mass nif, the pion mass m.,r, the vacuum expectation value of the a field a, the action S. 





Wj 


'"tt 


a 


S 


1.00 


.102(2) 


.22(1) 


.094(2) 


.363(2 


.95 


.133(2) 


.22(3) 


.128(2) 


.390(2 


.90 


.182(2) 


.19(3) 


.176(3) 


.428(2 


.85 


.233(3) 


.19(3) 


.231(3) 


.475(3 


.80 


.295(3) 


.21(3) 


.294(3) 


.535(3 


.75 


.359(2) 


.18(1) 


.365(3) 


.613(3 


.70 


.430(2) 


.18(2) 


.438(4) 


.707(3 


.65 


.503(3) 


.18(2) 


.519(4) 


.824(4 


.60 


.573(3) 


.18(1) 


.604(4) 


.969(5 


.55 






.695(4) 


1.14(1 


.50 


.746(2) 


.18(1) 


.795(3) 


1.36(1 
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TABLE II. A^/iatt = 3 data on 16^ lattice at coupling Xjg^ = 0.50 for various chemical 
potentials fi. Same notation as Table I but p is the fermion density and e is the energy density. 





ruf 




P 


(7 


S 


e 


.50 


.236(1) 


.188(2) 


.00045(2) 


.796(5) 


1.36(1) 


.217(1) 


.60 


.126(1) 


.24(6) 


.0069(4) 


.783(4) 


1.34(1) 


.223(1) 


.65 


.060(4) 


.21(1) 


.0163(4) 


.765(4) 


1.31(1) 


.232(1) 


.70 





.21(1) 


.0905(4) 


.586(4) 


1.10(2) 


.303(1) 


.725 




.40(2) 


.230(1) 


.131(4) 


.79(2) 


.426(1) 


.75 




.49(4) 


.267(1) 


.068(3) 


.77(1) 


.452(1) 


.775 




.55(8) 


.299(2) 


.050(3) 


.75(1) 


.475(1) 


.80 




.78(5) 


.337(2) 


.038(3) 


.74(1) 


.501(2) 


.825 




.78(5) 


.374(2) 


.026(2) 


.74(1) 


.526(2) 


.85 




— 


.416(2) 


.020(2) 


.73(1) 


.556(2) 






TABLE III. 


Same as Table II except = 0.80. 








ruf 




P 


a 


S 


e 


.15 


.137(2) 


.22(1) 


.0026(2) 


.287(1) 


.532(1) 


.303(1) 


.20 


.075(3) 


.20(1) 


.0047(1) 


.274(1) 


.525(1) 


.305(1) 


.225 


.041(4) 


.26(2) 


.0078(1) 


.262(2) 


.518(1) 


.308(1) 


.25 


.010(5) 


.22(1) 


.0120(7) 


.245(1) 


.511(1) 


.310(1) 


.275 





.24(1) 


.0162(6) 


.220(1) 


.501(1) 


.314(1) 


.30 




.25(2) 


.0220(3) 


.190(1) 


.490(1) 


.319(1) 


.325 




.28(3) 


.0288(2) 


.160(1) 


.480(1) 


.323(1) 


.35 




.30(3) 


.0375(4) 


.129(1) 


.472(1) 


.329(1) 


.40 




.35(3) 


.056(1) 


.083(1) 


.461(1) 


.339(1) 


.50 




.65(5) 


.099(1) 


.042(1) 


.452(1) 


.360(1) 
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TABLE IV. Same as Tahlo I except A'/iaii = 1. 





m/ 








a 


S 


e 


1.00 


.070(1) 




.25(5) 




.064(2) 


1.063(1) 


.314(1) 


.95 


.075(5) 




.22(3) 




.078(2) 


1.121(1) 


.313(1) 


.90 


.105(2) 




.23(3) 




.106(5) 


1.193(1) 


.309(1) 


.85 


.14(2) 




.20(2) 




.141(5) 


1.277(1) 


.305(1) 


.80 


1 no / /I ^ 
.19^(4) 




.21(2) 




.zOU(7j 


1.381(zj 


.z98(zj 


.70 


.32(1) 




.15(2) 




.336(9) 


1.655(2) 


.280(2) 


.60 


.475(5) 




.26(5) 




.511(12) 


2.065(2) 


.253(3) 


.oU 


.DOUl^O ) 




.15(1) 




.706(15) 


2.675(2) 


.zZ\j\o ) 




TART TT V 


S 3.1116 


as Table II except Nfi^tt = 1 and l/g"^ = 0.60. 






ruf 


m^r 




a 


P 


S 


e 


.30 


.173(1) 


.15(2) 




.498(2) 


.00133(57) 


2.062(1) 


.255(1) 


.35 


.116(2) 


.20(5) 




.491(2) 


.00432(19) 


2.055(1) 


.257(1) 


.40 


.046(2) 


.15(2) 




.470(2) 


.0104(3) 


2.038(1) 


.252(1) 


.425 









.439(4) 


.0171(4) 


2.018(2) 


.268(1) 


.45 




.19(2) 




.373(15) 


.0312(3) 


1.977(4) 


.281(1) 


.475 




.22(2) 




.259(15) 


.0541(8) 


1.924(3) 


.300(1) 


.4875 




.25(3) 




.206(15) 


.0644(3) 


1.904(2) 


.308(1) 


.50 




.30(2) 




.178(5) 


.0731(5) 


1.890(2) 


.313(1) 


.525 




.30(2) 




.124(3) 


.0892(4) 


1.871(1) 


.325(1) 


.55 




.36(3) 




.091(5) 


.1034(5) 


1.859(1) 


.333(1) 


.60 




.52(4) 




.060(7) 


.1310(7) 


1.843(1) 


.350(1) 
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TABLE VI. Same as Table V except l/g^ = 0.070. 





mj 




a 


P 


S 


e 


.15 


.164(1) 


.22(1) 


.331(1) 


.00067(33) 


1.65(1) 


.281(1) 


.20 


.113(5) 


.22(1) 


.324(1) 


.00247(64) 


1.65(1) 


.282(1) 


.225 


.085(5) 


.21(1) 


.315(1) 


.00468(9) 


1.64(1) 


.283(1) 


.25 


.040(8) 


.22(2) 


.301(1) 


.00689(43) 


1.63(1) 


.284(1) 


.275 




.24(1) 


.273(3) 


.0102(3) 


1.62(1) 


.289(1) 


.30 





.26(2) 


.235(7) 


.0175(4) 


1.61(1) 


.294(1) 


.325 




.30(1) 


.195(5) 


.0235(7) 


1.60(1) 


.299(1) 


.35 




.25(2) 


.162(7) 


.0309(6) 


1.59(1) 


.304(1) 


.40 




.38(1) 


.107(7) 


.0490(3) 


1.58(1) 


.314(1) 


.50 




.50(5) 


.052(3) 


.0898(3) 


1.56(1) 


.335(1) 
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Figure Captions 

Figure 1: Plot of (a) versus l/g"^ at fi = evaluated using the gap equation (2.14) on a 16^ 
lattice for three different bare fermion masses m. 

Figure 2: Gap equation prediction of (a) versus at l/g'^ = 0.5 for three different m. 
Figure 3: Gap equation prediction of (cr) versus fi ior m — 0.01 at four different 
Figure 4: Gap equation prediction of number density p versus p at three different plus 
the continuum free-field prediction. 

Figure 5: Gap equation prediction of energy density e versus /j, at three different 

Figure 6: Simulation results for a versus 1/g^ at p = using A^/iatt = 3 and m = 0.01. Also 

shown is the gap equation prediction of Fig. 1. 

Figure 7: Simulation results for a versus /j, a,t l/g"^ — 0.5 using A^/iatt = 3. Also shown is 
the gap equation prediction of Fig. 3. 

Figure 8: Simulation results for a versus m a,t l/g"^ — 0.5 and fj, — 0.5 using A^/iatt = 3. 
Figure 9: Simulation results for p and e versus p at 1/ g^ — 0.5 using A^/iatt — 3. 
Figure 10: Simulation results for action S versus p at l/g"^ = 0.5 using iVjiatt = 3. 
Figure 11: Simulation results for fermion mass m/ and pion mass at 1/g'^ — 0.5 using 
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Fig. 


11 except l/g"^ = 


= 0.8. 




Figure 


16: 


Same 


as 


Fig. 


6 except A^/iatt = 
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as 


Fig. 
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as 


Fig. 
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= 0.6. 
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as 


Fig. 
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= 1 and 1/g'^ 


= 0.6. 
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as 


Fig. 


11 except A^/iatt 
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= 0.6. 
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as 


Fig. 


17 except l/g^ = 


= 0.7. 




Figure 


22: 


Same 


as 


Fig. 


18 except l/g"^ - 


= 0.7. 
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Figure 23: Same as Fig. 19 except — 0.7. 
Figure 24: Same as Fig. 20 except l/g^ — 0.7. 
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